/*
NLA testsuite
(also contains some extra functions not reported in paper)

Arithmetic Algorithms
http://www.lsi.upc.edu/~erodri/webpage/polynomial_invariants/list.html : Enric Carbonell 'sinvariants
*/


void divbin_dummy(int A, int B,int q,int b, int r){}
int divbin (int A, int B){
  //binary division algorithm, by Kaldewaij
  //returns A//B
  //tvn: 2 whiles

  assert(A > 0);
  assert(B > 0);
  int q = 0;
  int r = A;
  int b = B;

  //inv(exists(int k;k >= 0  &&  b == 2^k*B)   &&   q == 0   &&   r == A  &&  A >= 0  &&  B > 0  );
  while(1){
    assert(q == 0);
    assert(r == A);
    assert(A >= 0);

    assert(B > 0);
    if (!(r>=b)) break;

    b=2*b;
    
  }

  //inv(A == q*b + r   &&   r >= 0   &&   exists(  int   k;   k >= 0 && b == 2^k*B  )   &&  B > 0  &&  b > r  );
  printf("A B q b r\n");
  while(1){
    assert(A == q*b + r);
    assert(r >= 0);
    assert(B > 0);
    assert(b > r);

    printf("%d %d %d %d %d\n",A,B,q,b,r);
    divbin_dummy(A,B,q,b,r);

    if (!(b!=B)) break;

    q = 2*q;

    b = b/2;
    if (r >= b){
      q = q + 1;
      r = r - b;
    }
  }
  
  assert(q==A/B);
  assert(r==A%B);
  
  return q;
}

void cohendiv_dummy(int x, int y, int q, int a, int b,int r){}
int cohendiv(int x, int y){
  //Cohen's integer division
  //returns x % y
  assert(x>0 && y>0);
  int q=0;
  int r=x;


  printf("x y q a b r\n");
  while(1){
    assert(x==q*y+r);
    assert(r>=0);
    if (!(r>=y)) break;

    int a=1;
    int b=y;
    while(1){
      assert(r>=y*a);
      assert(b==y*a);
      assert(x==q*y+r);
      assert(r>=0);
      printf("%d %d %d %d %d %d\n",x,y,q,a,b,r);
      cohendiv_dummy(x,y,q,a,b,r);

      if (!(r >= 2*b)) break;

      a = 2*a;
      b = 2*b;

    }
    r=r-b;
    q=q+a;
  }

  assert(r == x % y);
  assert(q == x / y);
  return q;
}


void mannadiv_dummy(y1,y2,y3,x1,x2){}
int mannadiv (int x1, int x2){
  //Division algorithm from
  //"Z. Manna, Mathematical Theory of Computation, McGraw-Hill, 1974"
  //return x1 // x2

  int y1,y2,y3;
  y1 = 0;
  y2 = 0;
  y3 = x1;
  

  printf("y1 y2 y3 x1 x2 \n");
  while(1){
    assert(y1* x2 + y2 + y3 == x1);
    printf("%d %d %d %d %d\n",y1,y2,y3,x1,x2);
    mannadiv_dummy(y1,y2,y3,x1,x2);
    
    if(!(y3 != 0))break;

    if (y2 + 1 == x2){
      y1 = y1 + 1;
      y2 = 0;
      y3 = y3 - 1;
    }
    else {
      y2 = y2 + 1;
      y3 = y3 - 1;
    }
  }

  assert(y1 == x1 / x2);
  return y1;
}


void hard_dummy(int d1, int A, int D, int p, int r){}
int hard(int A, int D){
  //hardware integer division program, by Manna
  //returns q==A//D

  int r,d1,p,q;

  r=A;
  d1=D;
  p=1;
  q=0;
  
  //inv( A >= 0 && D > 0 && q == 0 && r == A && d1 == D*p && exists( int k; k>=0 && d1==D*2^k) && exists( int l; l>=0 && p==2^l ) );
  while(1){
    assert(A >= 0);
    assert(D > 0);
    assert(q == 0);
    assert(r == A);
    assert(d1 == D*p);

    if(!(r>= d1)) break;
    d1=2*d1;
    p=2*p;
  }

  //inv( D > 0 && A == q*D+r && d1 > r && r >= 0 && d1 == D*p && exists( int k; k>=0 && d1==D*2^k) && exists( int l; l>=0 && p==2^l ) );
  printf("d1 A D q p r\n");
  while(1){
    assert(D > 0);
    assert(A == q*D+r);
    assert(d1 > r);
    assert(r >= 0);
    assert(d1 == D*p);
    printf("%d %d %d %d %d %d\n", d1, A,D, q, p, r);
    hard_dummy(d1, A, D, p, r);

    if(!(p!=1))break;

    d1=d1/2;
    p=p/2;
      
    if(r>=d1){
      r=r-d1;
      q=q+p;
    }

  }

  //post( N == q*D+r && r >= 0 && r < D );
  assert(A = q*D+r);
  assert(r >= 0);
  assert(r < D);
  assert(r == A % D);
  assert(q == A / D);

  return q;
}




void wensley0_dummy(float P, float Q, float E, 
		    float a, float b, float d, float y){}
float wensley0 (float P, float Q, float E){
  /* Wensley's algorithm for real division from Rodriguez's website*/
  //pre(Q > P  &&  P >= 0  &&  E > 0   );
  assert(Q > P);
  assert(P >= 0);
  assert(E > 0);

  float a,b,d,y;
  a=0;
  b=Q/2;
  d=1;
  y=0;
  
  //inv(a == Q*y   &&   b == Q*d/2   &&   P/Q - d < y   &&   y <= P/Q   );
  printf("P Q E a b d y\n");
  while(1){
    assert(a == Q*y);
    assert(b == Q*d/2);
    assert(P/Q - d < y);
    assert(y <= P/Q);
    printf("%f %f %f %f %f %f %f\n",P,Q,E,a,b,d,y);
    wensley0_dummy(P,Q,E,a,b,d,y);

    if(!(d>= E))break;

    if (P< a+b){
      b=b/2;
      d=d/2;
      }
    else{
      a=a+b;
      y=y+d/2;
      b=b/2;
      d=d/2;
    }
  }
  //post(P/Q >= result   &&   result > P/Q - E   );
  assert(P/Q >= y);
  assert(y > P/Q - E);
  return y;

}


//tvn: todo :  rerurn wensley with new set of inputs
void wensley_dummy(float x, float y, float t, float z, 
                   float d, float u, float v){}
float wensley (float x, float y, float t){
  /* Wensley's algorithm for real division */
  //from http://www.cs.indiana.edu/classes/c241-john/post/Wensley.pdf
  /*
    some inputs don't work properly due to floating point imprecision
    e.g., wensley 0.85 0.93 0.02
    u 0.813750, z 0.875000, y 0.930000
    even though  0.813750 == 0.875000 *  0.93000,
    computer gives 0.813750 < 0.875000 *  0.93000
   */
  
  assert(0 <= x);
  assert(x<y);
  assert(y <= 1);
  
  float z,d,u,v;
  z = 0.0;
  d = 1.0;
  u = 0.0;
  v = y/2.0;


  printf("x y t z d u v\n");
  while(d > t){
    //yz <= x
    //u == z*y  #should be == instead of  u>= z*y
    //v == d/2.0*y

    printf("%f %f %f %f %f %f %f\n",x, y, t, z, d, u, v);
    wensley_dummy(x, y, t, z, d, u, v);

    d = d/2.0;

    if (u + v <= x){
      z = z + d; 
      u = u + v; 
    }
    v = v/2.0;
    //2v == d*y
  }

  assert(z <= x/y &&  x/y < z + t);
  return z;
}//wensley



void sqrt1_dummy(int a, int n,int t,int s){}
int sqrt1(int n){
  /* computing the floor of the square root of a natural number */
  assert(n>=0);

  int a,s,t;
  a=0;
  s=1;
  t=1;

  printf("a n t s\n");
  //inv( a ^2 <= n   &&   t == 2*a + 1   &&   su == ( a + 1 ) ^2    );    
  while(1){

    assert((int)pow(a,2) <= n);
    assert(t == 2*a + 1);
    assert(s == (int)pow((a + 1),2));
    printf("%d %d %d %d\n",a,n,t,s);
    sqrt1_dummy(a,n,t,s);

    if(!( s <= n )) break;

    a=a+1;
    t=t+2;
    s=s+t;

  }
  //post(result ^2 <= n    &&( result + 1 ) ^2 > n   ); 
  assert((int)pow(a,2) <= n);
  assert((int)pow((a+1),2) > n);
  assert(a == (int)floor(sqrt(n)));
  return a;
}//sqrt1

void freire1_dummy(float a,float x,int r){}
int freire1 (float a ){
  //algorithm for finding the closest integer to the square root,
  //from  www.pedrofreire.com/crea2_en.htm?

  float x;
  int r;

  x=a/2.0;
  r=0.0;

  
  printf("a x r\n");
  //inv( a == 2*x + r^2 - r   &&   x>0 );

  while(1){
    assert(a == 2*x + (int)pow(r,2) - r);
    assert( x>0);
    printf("%f %f %d\n",a,x,r);
    freire1_dummy(a,x,r);

    if(!( x>r )) break;

    x = x-r;
    r = r+1.0;
    
  }

  //  post( result^2 + result >= a  &&  result^2 - result < a );
  assert((int)pow(r,2)+r >= a);
  assert((int)pow(r,2)-r < a);
  assert(r == (int)round(sqrt(a)));

  return r;
}//freire1

void freire2_dummy(float a,float x,float s, int r){}
int freire2(float a){
  //algorithm for finding the closest integer to the cubic root,
  //from www.pedrofreire.com/crea2_en.htm? */ 
  //a is a float

  assert(a>0);

  float x,s;
  int r;
  
  x=a;
  r=1;
  s=3.25;

           
  printf("a x s r\n");
  //inv(4*r^3 - 6*r^2 + 3*r + 4*x - 4*a == 1   &&   x>0   &&   -12*r^2 + 4*s== 1   );

  while(1){
    assert(4*(int)pow(r,3) - 6*(int)pow(r,2) + 3*r + 4*x - 4*a == 1);
    assert(x>0);
    assert(-12*(int)pow(r,2) + 4*s== 1);

    printf("%f %f %f %d\n",a,x,s,r);
    freire2_dummy(a,x,s,r);
    if(!(x-s > 0.0)) break;

    x = x - s;
    s = s + 6 * r + 3;
    r = r + 1;
  }

  //post( 4*r^3+6*r^2+3*r >= 4*a    &&   4*a > 4*r^3-6*r^2+3*r-1     );
  assert(4*(int)pow(r,3)+6*(int)pow(r,2)+3*r >= 4*a);
  assert(4*a > 4*(int)pow(r,3)-6*(int)pow(r,2)+3*r-1);
  assert(r == (int)round(pow(a,(1.0/3.0))));
  return r;
}//freire2

void dijkstra_dummy(int r, int p, int n, int q, int h){}
int dijkstra(int n){
  /* program computing the floor of the square root, by Dijkstra */
  int p,q,r,h;

  

  p=0;
  q=1;
  r=n;
  h = -12345;//dummy initial value
  
  //inv(n >= 0   &&   p == 0   &&   r==n   &&  exists ( int   k; k >= 0  &&  q == 4^k )  );
  while(1){
    assert(n >= 0);
    assert(p == 0);
    assert(r == n);

    if(!(q<=n))break;

    q=4*q;
  }

  //inv(  exists (  int  l;  l >= 0  &&  q ==4^l )   &&   r >= 0   &&   r < 2*p + q   &&   p^2 + r*q == n*q   ); 
  //printf("#inv: r < 2*p + q,  p*p + r*q == n*q\n");
  printf("r p n q h\n");
  while(1){
    assert(r >= 0);
    assert(r < 2*p + q);
    assert((int)pow(p,2) + r*q == n*q);

    printf("%d %d %d %d %d\n",r,p,n,q,h);
    dijkstra_dummy(r,p,n,q,h);

    if (!(q!=1)) break;

    q=q/4;
    h=p+q;
    p=p/2;
    if (r>=h){
      p=p+q;
      r=r-h;
    } 

  }
  //post(  result^2 <= n   &&( result + 1) ^2 > n  );
  assert((int)pow(p,2)<=n);
  assert((int)pow(p+1,2)>n);
  assert(p == (int)floor(sqrt(n)));
  return p;
}


void cohencu_dummy(int a, int n, int x,int y,int z){}
int cohencu(int a){
  /* printing consecutive cube, by Cohen */

  //pre( true );

  int n,x,y,z;
  int cubes[a];

  n=0; x=0; y=1; z=6;
  
  printf("a n x y z\n");
  //inv(  z == 6*n + 6   &&   y == 3*n^2 + 3*n + 1   &&   x == n^3  );
  while(1){
    assert(z == 6*n + 6);
    assert(y == 3*(int)pow(n,2) + 3*n + 1);
    assert(x == (int)pow(n,3));

    printf("%d %d %d %d %d\n", a,n,x,y,z);
    cohencu_dummy(a,n,x,y,z);

    if(!(n<=a)) break;
    
    assert(x==(int)pow(n,3));
    cubes[n]=x;
    n=n+1;
    x=x+y;
    y=y+z;
    z=z+6;
  }

  print_array(cubes,a,0);
  return -1;
}//cohencu


void euclidex1_dummy(int a , int b , int c ,int p , 
                     int q, int r ,int s, int x ,int y ,
                     int k){}    

int euclidex1 (int x, int y){
  /* extended Euclid's algorithm */
  assert(x > 0);
  assert(y > 0); 
  int a,b,p,q,r,s;


  a=x;
  b=y;
  p=1;
  q=0;
  r=0; 
  s=1;


  printf("a b c p q r s x y k\n");
  //inv( GCD(a,b)==GCD(x,y) && a==y*r+x*p && b==x*q+y*s );
  while(1){ 
    assert(a == y*r+x*p);
    assert(b == x*q+y*s);

    if(!(b!=0)) break;
    int c,k;
    c=a;
    k=0;

    //inv( a == k*b+c && GCD(a,b) == GCD(x,y) && a==y*r+x*p && b==x*q+y*s );
    while(1){
 
      assert(a == k*b+c);
      assert(a == y*r+x*p);
      assert(b == x*q+y*s);

      printf("%d %d %d %d %d %d %d %d %d %d\n",
             a ,b ,c ,p , q,r ,s, x ,y ,k);
      euclidex1_dummy(a ,b ,c ,p , q,r ,s, x ,y ,k);

      if(!( c >= b )) break;

      c=c-b;
      k=k+1;

    }


    a=b;
    b=c;

    int temp;
    temp=p;
    p=q;
    q=temp-q*k;
    temp=r;
    r=s;
    s=temp-s*k;
  }
  assert(mygcd(x,y)==a);
  assert(a == y*r+x*p);
  return a;
}//euclidex1


void euclidex2_dummy(int a, int b,int y,int r,int x,int p,int q, int s){}
int euclidex2 (int A, int B){
  /* extended Euclid's algorithm */
  assert(A > 0);
  assert(B > 0); 

  int x,y,i,j,k,m;
    
  x=A;
  y=B;
  i=1;
  j=0;
  k=0;
  m=1;
  
  printf("A B x y i j k m\n");
  //inv( GCD(x,y) == GCD(A,B)    &&    x == B*k + A*i    &&    y == A*j + B*m    );
  while(1){ 
    assert(mygcd(x,y) == mygcd(A,B));
    assert(x == A*i + B*j);
    assert(y == A*k + B*m);

    printf("%d %d %d %d %d %d %d %d\n",A,B,x,y,i,j,k,m);
    euclidex2_dummy(A,B,x,y,i,j,k,m);

    if(!(x!=y)) break;

    if (x>y){x = x-y; i = i-k; j=j-m;}
    else {y = y-x; k = k-i; m = m-j;}

  }
  assert(mygcd(A,B) == x);
  assert(x == A*i + B*j   );

  return x;
}//euclidex2


void euclidex3_dummy(int a,int b,int y,int r,int x,int p,
                     int q,int s,int d,int v, int k, int c){}
int euclidex3 (int x, int y){
  /* extended Euclid's algorithm */
  assert(x > 0);
  assert(y > 0); 

  int a,b,p,q,r,s;

  a=x; b=y;  p=1;  q=0;  r=0;   s=1;

  assert(a==y*r+x*p);
  assert(b==x*q+y*s);


  printf("a b y r x p q s d v k c\n");
  //inv( GCD(a,b)==GCD(x,y) && a==y*r+x*p && b==x*q+y*s );
  while(1){ 
    assert(mygcd(a,b) == mygcd(x,y));
    assert(a == y*r+x*p);
    assert(b == x*q+y*s);

    if(!(b!=0)) break;

    int c,k;
    c=a;
    k=0;

    //inv( GCD(a,b)==GCD(x,y) && a==y*r+x*p && b==x*q+y*s && a==quot*b+c ); 
    while(1){
      assert (mygcd(a,b) == mygcd(x,y));
      assert(a == y*r+x*p);
      assert(b == x*q+y*s);
      assert(a == k*b+c ); 
      
      if(!(c>=b)) break;

      int d,v;
      d=1;
      v=b;
    
      //inv( GCD(a,b)==GCD(x,y) && a==y*r+x*p && b==x*q+y*s && a==k*b+c && v==b*d ); 
      while(1){
	assert(mygcd(a,b) == mygcd(x,y));
	assert(a == y*r+x*p);
	assert(b == x*q+y*s);
	assert(a == k*b+c);
	assert(v == b*d);
	
        printf("%d %d %d %d %d %d %d %d %d %d %d %d\n", 
               a, b, y, r, x, p, q, s, d, v, k, c);
        euclidex3_dummy(a, b, y, r, x, p, q, s, d, v, k, c);

	if(!(c>= 2*v)) break;

        d = 2*d;
        v = 2*v;

      }
      c=c-v;
      k=k+d;

      //assert(a==y*r+x*p && b==x*q+y*s && a==k*b+c );
             
    }
      
    a=b;
    b=c;
    int temp;
    temp=p;
    p=q;
    q=temp-q*k;
    temp=r;
    r=s;
    s=temp-s*k;
  }
      
  assert(mygcd(x,y)==a);
  assert(a==y*r+x*p );
      
  return a;
}//euclidex3
    

void lcm1_dummy(int a, int b, int x, int y, int u, int v){}
int lcm1(int a, int b){
  /* algorithm for computing simultaneously the GCD and the LCM, 
     by Sankaranarayanan */
  /*tvn: 2 whiles*/
  assert(a>0);
  assert(b>0);

  int x,y,u,v;

  x=a;
  y=b;
  u=b;
  v=0;
  
  //inv( GCD(x,y) == GCD(a,b) && x*u + y*v == a*b );
  printf("a b x y u v\n");
  while(1){

    //inv(GCD(x,y) == GCD(a,b) && x*u + y*v == a*b );
    assert(mygcd(x,y) == mygcd(a,b));
    assert (x*u + y*v == a*b);
    printf("%d %d %d %d %d %d\n", a,b,x,y,u,v);
    lcm1_dummy(a,b,x,y,u,v);
    
    if(!(x!=y)) break;

    //inv( GCD(x,y) == GCD(a,b) && x*u + y*v == a*b ); 
    while(1){
      assert(mygcd(x,y) == mygcd(a,b));
      assert(x*u + y*v == a*b); 
      if(!(x>y)) break;
      x=x-y;
      v=v+u;
    }
    
    //inv( GCD(x,y) == GCD(a,b) && x*u + y*v == a*b );
    while(1){
      assert(mygcd(x,y) == mygcd(a,b));
      assert(x*u + y*v == a*b);
      if(!(x<y)) break;
      y=y-x;
      u=u+v;
    }

  }
  assert(x == mygcd(a,b));
  assert(u+v==mylcm(a,b));
  return u+v; //lcm
}


void lcm2_dummy(int a,int b,int x,int y,int u,int v){}
int lcm2 (int a, int b){
  /* algorithm for computing simultaneously the GCD and the LCM, by Dijkstra */
  assert(a>0);
  assert(b>0);

  int x,y,u,v;

  x=a;
  y=b;
  u=b;
  v=a;

  printf("a b x y u v\n");
  //inv(GCD(x,y) == GCD(a,b)   &&   x>0   &&   y>0   &&   x*u + y*v == 2*a*b  );
  while(1){ 
    assert(mygcd(x,y) == mygcd(a,b));
    assert(x>0);
    assert(y>0);
    assert(x*u + y*v == 2*a*b);

    printf("%d %d %d %d %d %d\n", a,b,x,y,u,v);
    lcm2_dummy(a,b,x,y,u,v);

    if(!(x!=y)) break;
    if (x>y){
      x=x-y;
      v=v+u;
    }
    else {
      y=y-x;
      u=u+v;
    }
    
  }
  assert(x == mygcd(a,b));
  assert((u+v)/2 == mylcm(a,b));
  return (u+v)/2; //lcm
}//lcm2


void prodbin_dummy(int a, int b, int x, int y, int z){}
int prodbin (int a,int b){
  /* algorithm for computing the product of two natural numbers */
  /* shift_add */

  assert(a>=0);
  assert(b>=0);

  int x,y,z;
  
  x = a;
  y = b;
  z = 0;
  
  printf("a b x y z\n");
  //inv(z+x*y==a*b );
  while(1){ 
    assert(z+x*y==a*b);
    printf("%d %d %d %d %d\n", a,b,x,y,z);
    prodbin_dummy(a,b,x,y,z);

    if(!( y!=0 )) break;

    if (y%2 ==1 ){
      z = z+x;
      y = y-1;
    }
    x = 2*x;
    y = y/2;


  }
  assert(z == a*b);
  return z; 
}//prodbin


void prod4br_dummy(int x, int y, int a, int b, int p, int q){}
int prod4br(int x, int y){
  /* algorithm for computing the product of two natural numbers */
  int a,b,p,q;
  assert(x>=0);
  assert(y>=0);
  
  a = x;
  b = y;
  p = 1;
  q = 0;

  printf("x y a b p q\n");
  //inv( q+abp==xy );
  while(1){ 

    assert(q+a*b*p==x*y);
    printf("%d %d %d %d %d %d\n", x, y, a, b, p, q);
    prod4br_dummy(x, y, a, b, p, q);

    if(!(a!=0 && b!=0 )) break;

    if (a % 2 ==0 && b % 2 ==0 ){
      a =a/2;
      b = b/2;
      p = 4*p;
    }
    else if (a % 2 ==1 && b % 2 ==0 ){
      a =a-1;
      q = q+b*p;
    }
    else if (a % 2 ==0 && b % 2 ==1 ){
      b =b-1;
      q = q+a*p;
    }
    else {
      a =a-1;
      b =b-1;
      q = q+(a+b+1)*p;  /*dammit I am good ---  was (a+b-1)*/
    }
  }

  assert(q == x*y);
  return q; 
}

void fermat1_dummy(int A, int R, int u, int v, int r){}
int fermat1(int A, int R){
  // program computing a divisor for factorisation, by Knuth 4.5.4 Alg C ?
  /*tvn: run forever, but ok because we can just capture those values*/
  //tvn: 2 whiles

  //pre( A>=1 && (R-1)*(R-1) < A && A<=R*R && A%2==1 );
  assert(A>=1);
  assert((R-1)*(R-1) < A);
  assert(A<=R*R);
  assert(A%2==1);

  int u,v,r;
  u=2*R+1;
  v=1;
  r=R*R-A;
  
  //inv( 4*(A+r)==u*u-v*v-2*u+2*v && v%2==1 && u%2==1 && A>=1 );

  printf("A R u v r\n");
  while(1){
    
    assert(4*(A+r)==u*u-v*v-2*u+2*v);
    assert(v%2==1);
    assert(u%2==1);
    assert(A>=1);

    printf("%d %d %d %d %d\n", A, R, u, v, r);
    fermat1_dummy(A, R, u, v, r );    

    if(!( r!=0)) break;


    //inv( 4*(A+r)==u*u-v*v-2*u+2*v && v%2==1 && u%2==1 && A>=1 ); 
    while(1){
      assert(4*(A+r)==u*u-v*v-2*u+2*v);
      assert(v%2==1);
      assert(u%2==1);
      assert(A>=1); 

      if(!(r>0)) break;

      r=r-v;
      v=v+2;
    }
    
    //inv( 4*(A+r)==u*u-v*v-2*u+2*v && v%2==1 && u%2==1 && A>=1 );  
    while(1){
      assert(4*(A+r)==u*u-v*v-2*u+2*v);
      assert(v%2==1);
      assert(u%2==1);
      assert(A>=1);

      if(!(r<0)) break;
      r=r+u;
      u=u+2;
    }

  }
  
  assert(u!=v );
  assert(A % ((u-v)/2) == 0);
  return((u-v)/2);

}//fermat1

void fermat2_dummy(int A, int R, int u, int v, int r){}
int fermat2(int A, int R){
  /* program computing a divisor for factorisation, by Bressoud */

  //pre(  A >= 1  &&(R-1) ^2 < A   &&   A <= R^2  &&  A%2 == 1  );
  assert(A >= 1);
  assert((R-1)*(R-1) < A);
  assert(A<=R*R);
  assert(A%2 == 1);

  int u,v,r;
    
  u=2*R+1;
  v=1;
  r=R*R-A;


  printf("A R u v r\n");
  //inv(4*(A + r) == u^2 - v^2 - 2*u + 2*v   &&   v%2 == 1   &&   u%2 == 1   &&   A >= 1  );
  while(1){

    assert(4*(A + r) == u*u - v*v - 2*u + 2*v);
    assert(v%2 == 1);
    assert(u%2 == 1);
    assert(A >= 1);

    printf("%d %d %d %d %d\n", A, R, u, v, r);
    fermat2_dummy(A, R, u, v, r);
    
    if(!(r!=0)) break;
    
    if (r>0){
      r=r-v;
      v=v+2;
    }
    else{
      r=r+u;
      u=u+2;
    }
  }
  
  assert(u!=v);
  assert(A % ((u-v)/2) == 0);
  return((u-v)/2);

}//fermat2

void knuth_dummy(int n, int d1, int r, int k, int q, int d, 
                 int s, int t){}
int knuth(int n, int d1){
  //algorithm searching for a divisor for factorization, by Knuth
  //had to comment out one of the invariants
  //pre( n % 2 == 1    &&    n > 0    &&    d1 % 2 == 1    && ( d1 - 1 ) ^3 < 8*n    );
  assert(n % 2 == 1);
  assert(n > 0);
  assert(d1 % 2 == 1);
  assert((int)pow((d1 - 1),3) < 8*n);

  int r,k,q,d,s,t;

  d=d1;
  r= n % d;

  if (d <= 2){
    printf("#E: d - 2 <= 0.  Use a diff val for d\n");
    return 0;
  }
  k=n % (d-2);
  q=4*(n/(d-2) - n/d);
  s=sqrt(n);
  

  printf("n d1 r k q d s t\n");
  //inv(d^2*q - 4*r*d + 4*k*d - 2*q*d + 8*r == 8*n    &&    d % 2 == 1   );

  while(1){

    assert(d*d*q - 4*r*d + 4*k*d - 2*q*d + 8*r == 8*n);
    assert(d % 2 == 1);

    printf("%d %d %d %d %d %d %d %d\n", n, d1, r, k, q, d, s, t);
    knuth_dummy(n, d1, r, k, q, d, s, t);

    if(!((s>=d)&&(r!=0))) break;

    if (2*r-k+q<0){
      t=r;
      r=2*r-k+q+d+2;
      k=t;
      q=q+4;
      d=d+2;
    } 
    else if ((2*r-k+q>=0)&&(2*r-k+q<d+2)){
      t=r;
      r=2*r-k+q;
      k=t;
      d=d+2;
    }
    else if ((2*r-k+q>=0)&&(2*r-k+q>=d+2)&&(2*r-k+q<2*d+4)){
      t=r;
      r=2*r-k+q-d-2;
      k=t;
      q=q-4;
      d=d+2;
    }
    else {/* ((2*r-k+q>=0)&&(2*r-k+q>=2*d+4)) */
      t=r;
      r=2*r-k+q-2*d-4;
      k=t;
      q=q-8;
      d=d+2;
    }

  }
  // assert(r == 0 => ( n % d == 0 )   );
  assert(!(r==0) ||  ( n % d == 0 )  );
  return d;

}//knuth




void petter(k){  
  //geo
  int x = 0;
  int y = 0;
  printf("x y k\n");
  while(1){
    printf("%d %d %d\n",x,y,k);
    if (!(y != k )) break;
    x = x + pow(y,5);
    y = y + 1;
    
  }
}


/*from petter"s master thesis
  http://www2.cs.tum.edu/~petter/da/da.pdf*/
int geo1 (int z, int k){
  //todo: probably not use, since no loop inv
  assert(k>0);
  int x = 1;
  int y = z;
  int k1 = 1;


  printf("x y z\n");
  while(1){
    //1 + x - y == 0
    assert(y == pow(z,k));    
    printf("%d %d %d\n",x,y,z);

    if(!(k1 < k)) break;

    k1 = k1 + 1;
    x = x*z + 1;
    y = z*y;

    x = x *(z -1);


  }//geo1



  //note the paper says x = (z-1)* sum(z^k)[k=0..k] is not correct
  //this code produces  x = (z-1)* sum(z^k)[k=0..k-1]
  //assert(x == (z-1) * sum([z**k_ for k_ in range(k)]));
  //return x,y
  return x;
}

void geo2_dummy(int x, int y, int z, int k){}
int geo2(int z, int k){
  assert (k>0);
  int x = 1;
  int y = 1;
  int c = 1;


  printf("x y z k\n");
  while(1){
    // 1+x*z-x-z*y==0
    printf("%d %d %d %d\n",x, y, z, k);
    geo2_dummy(x, y, z, k);

    if(!(c < k)) break;

    c = c + 1;
    x = x*z + 1;
    y = y*z;
  }
  
  assert(y == pow(z,(k-1)));
  //assert(x == sum([z**k_ for k_ in range(k)]))
  //return x,y
  return x;

}//geo2

void geo3_dummy(int x, int y, int z, int a, int k){}
int geo3(int z,int a, int k){
  int x = a;
  int y = 1;
  int c = 1;


  printf("x y z a k\n");
  while(1){

    //z*x-x+a-a*z*y == 0
    printf("%d %d %d %d %d\n", x, y, z, a,k) ;
    geo3_dummy( x, y, z, a,k);
    if(!(c < k)) break;

    c = c + 1;
    x = x*z + a;
    y = y*z;
  }

  assert(y == pow(z,(k-1)));
  //assert(x == sum([a*(z**k_) for k_ in range(k)]));
  return x;
}//geo3


/*
  # // yielding x=y
  # //todo: too easy, not worth doing
  # def ps1 ( String []
  # y = 0
  # x = 0
  # while( args !=null){
  # y=y +1
  # x=x +1
*/

void ps2_dummy(int x, int y, int k){}
int ps2(int k){
  int y = 0;
  int x = 0;
  int c = 0;


  printf("x y k\n");
  while(1){
    //2*x-pow(y,2)-y == 0

    printf("%d %d %d\n",x,y,k);
    ps2_dummy(x,y,k);
    if(!(c < k))break;

    c = c + 1;
    y=y +1;
    x=x+y;

  }

  assert(y == 1*k);
  //assert(x == sum(range(k+1)))
  
  return x;
}//ps2

void ps3_dummy(int x, int y, int k){}
int ps3(int k){
  int y = 0;
  int x = 0;
  int c = 0;
  

  printf("x y k\n");
  
  while(1){
    //6*x-2*pow(y,3)-3*pow(y,2)-y == 0

    printf("%d %d %d\n",x,y,k);
    ps3_dummy(x,y,k);
    if(!(c < k))break;
    c = c + 1;
    y=y +1;
    x=y*y+x;

  }

  assert(y == 1*k);

  //assert(x == sum(k_**2 for k_ in range(k+1)))
  return x;
}//ps3
  

void ps4_dummy(int x,int y,int k){}
int ps4 (int k){
  int y = 0;
  int x = 0;
  int c = 0;
  
  printf("x y k\n");
  
  while(1){
    //4*x-y**4-2*y**3-y**2 == 0
    printf("%d %d %d\n",x,y,k);
    ps4_dummy(x,y,k);
    if(!(c < k))break;

    c = c +1 ;
    y=y +1;
    x=y*y*y+x;
  }

  assert(y == 1*k);
  //assert(x == sum(k_**3 for k_ in range(k+1)))
  //return x,y
  return x;
}//ps4



void readers_writers_dummy(int r, int w, int k, int a1, int a2, int k0){}
int readers_writers(int a1, int a2, int k0){
  //generalization of the readers-writers problem, by Sankaranarayanan
  assert ( a1 >= 0 && a2 >= 0 && k0 >=0 ) ;
  int r,w,k;
  r = 0;
  w = 0;
  k = k0;

  int ctr = 0;

  printf("r w k a1 a2 k0\n");
  while(ctr < MAX_LOOP){ 
    
    // r*a1 + w*a2 + k == k0
    //r*w==0
    printf("%d %d %d %d %d %d\n",r, w, k, a1, a2, k0);
    readers_writers_dummy(r, w, k, a1, a2, k0);

    ctr++;

    if(w==0){
      r = r+1;
      k = k-a1;
    }
    if(r==0){
      w = w+1;
      k = k-a2;
    }
    if( w==0){
      r = r-1;
      k = k+a1;
    }
    if(r==0){
      w = w-1;
      k = k+a2;
    }

  }
  return k;
}//readers_writers

void z3sqrt_dummy(float a, float err, float r, float p, float q){}
float z3sqrt (float a, float err){
  //program for computing square roots, by Zuse

  //todo: due to rounding imprec
  //assert(1.0 <= a && a <= 4.0 && 0.005 <= err );

  float r,q,p;

  r = a-1.0;
  q = 1.0;
  p = 1.0/2.0;
  
  printf("a err r p q\n");

  //inv( q^2 + 2*r*p == a && err > 0 && p >= 0 && r >= 0 );
  while(1){
    
    assert(q*q + 2*r*p == a);
    assert(err > 0);
    assert(p >= 0);
    assert(r >= 0 );

    printf("%f %f %f %f %f\n",a,err,r,p,q);
    z3sqrt_dummy(a,err,r,p,q);

    if(!(2.0*p*r >= err)) break;

    if ( 2.0*r-2*q-p >= 0.0 ){
      r = 2.0*r-2.0*q-p;
      q = q+p;
      p = p/2.0;
    }
    else{
      r = 2.0*r;
      p = p/2.0;
    }
    
    //tvn:  changed to 'q*q + 2.0*r*p >= a' from '=='
    //assert(q*q + 2.0*r*p == a);
    //assert(err > 0.0);
    //assert(p >= 0.0);
    //assert(r >= 0.0 );
    
  }

  assert(pow(q,2) <= a);
  assert(pow(q,2)+err > a);
  return q;
  
}//z3



int f1(int a){
  assert (a >= 0);
  int b = 0; 
  int c = 1;

  while(c*c <= a){
    c = 2*c;
  }

  while(c >= 2){
    
    c = c/2;
    if ((int)pow((b + c),2.0) <= a){
      assert ((int)pow((2*b + c),2.0) <= 4*a);
      b = b + c;
    }
    else{
      printf("%d %d %d, %d , %d\n",a,b,c,(int)pow((2*b + c),2.0),4*a+4);
      assert ((int)pow((2*b + c),2.0) >= 4*a+4);
    }
  }

  assert(b == (int)floor(sqrt(a)));
  return b;
}
